__global__ void Newton(unsigned char* surfaceIn, unsigned char* surfaceOut, size_t pitch, uint pixelBits, uint patchBits, 
								  uint* PointOccupiedSortPad, uint* ErrorOccupied, uint* ErrorOccupiedValid, uint* BoundaryOccupied, uint* BoundaryOccupiedValid,
								  float* FuncData, int* PatchIndex, int* PatchIndexInverse, int* PolarIndex, ABSSURF_SHARE_INFO asi, int iter)
{	
	int ind = threadIdx.x;
	__shared__ uint patchInfo[PIXELTHREADS];
	patchInfo[ind] = PointOccupiedSortPad[ind+blockIdx.x*blockDim.x];
	__syncthreads();

	int patchIndex = (int)( (patchInfo[0]<<(pixelBits+pixelBits))>>(pixelBits+pixelBits) );
	int segInd = PatchIndexInverse[patchIndex];
	//if(0!=patchIndex)	return;
	int patchIndexXYZ[3];
	patchIndexXYZ[0] = segInd/asi.segShift[0];
	patchIndexXYZ[1] = segInd%asi.segShift[0]/asi.segShift[1];
	patchIndexXYZ[2] = segInd%asi.segShift[1];
	float patchXYZ[3];
	patchXYZ[0] = patchIndexXYZ[0]*asi.boxSize[0] + asi.box0[0];
	patchXYZ[1] = patchIndexXYZ[1]*asi.boxSize[1] + asi.box0[1];
	patchXYZ[2] = patchIndexXYZ[2]*asi.boxSize[2] + asi.box0[2];

	bool useRegulaFalsi = false;
	if(PolarIndex[patchIndex]>0)	useRegulaFalsi = true;
	//patchIndex = PatchIndex[patchIndex];
    //if (patchIndex<0)                           return;


	extern __shared__ float funcNormal[];
	int funcSize = (asi.deg[0]+1)*(asi.deg[1]+1)*(asi.deg[2]+1);
	float* func = funcNormal;
#ifdef NORMALFUNCTION
	float* funcJacobi = &func[funcSize];
#else
#endif


	//if(ind<funcSize)
		func[ind] = FuncData[patchIndex*funcSize+ind];
	__syncthreads();


#ifdef NORMALFUNCTION
	int degShift[3] = {(asi.deg[1]+1)*(asi.deg[2]+1), (asi.deg[2]+1), 1};
	if((ind/degShift[0])<asi.deg[0])
	{
		funcJacobi[ind]				= asi.deg[0]*(func[ind +degShift[0]] - func[ind] );
	}
	if((ind%degShift[0]/degShift[1])<asi.deg[1])
	{
		funcJacobi[funcSize +ind]	= asi.deg[1]*(func[ind +degShift[1]] - func[ind] );
	}
	if((ind%degShift[1])<asi.deg[2])
	{
		funcJacobi[funcSize*2 +ind] = asi.deg[2]*(func[ind +degShift[2]] - func[ind] );
	}
	__syncthreads();
#else
#endif


	uint pixelX = (patchInfo[ind]>>(patchBits+pixelBits));
	uint pixelY = (patchInfo[ind]<<pixelBits)>>(patchBits+pixelBits);
	float* pixel = (float*)(surfaceOut + pixelY*pitch) + 4*pixelX;

	if(pixel[3]!=0.0f)			return;

	if(patchInfo[ind]==0)		return;


	float xyz[3] = {0.0f, 0.0f, 0.0f};
	float norm[3] = {0.0f, 0.0f, 0.0f};
	float fValue[3] = {0.0f, 0.0f, 0.0f};
	float nablaF[3] = {0.0f, 0.0f, 0.0f};
	float theta[3];


	float screenX = 2.0f*(pixelX+0.5f)/cRenderParam.renderWidth-1.0f;
	float screenY = 1.0f-2.0f*(pixelY+0.5f)/cRenderParam.renderHeight;
	float screenZ = *((float*)(surfaceIn+pixelY*pitch) + pixelX*4 + 2);

	float worldW = screenX*cRenderParam.invViewProj[0][3] + screenY*cRenderParam.invViewProj[1][3] + screenZ*cRenderParam.invViewProj[2][3] + cRenderParam.invViewProj[3][3];
	xyz[0] = (screenX*cRenderParam.invViewProj[0][0] + screenY*cRenderParam.invViewProj[1][0] + screenZ*cRenderParam.invViewProj[2][0] + cRenderParam.invViewProj[3][0]) / worldW;
	xyz[1] = (screenX*cRenderParam.invViewProj[0][1] + screenY*cRenderParam.invViewProj[1][1] + screenZ*cRenderParam.invViewProj[2][1] + cRenderParam.invViewProj[3][1]) / worldW;
	xyz[2] = (screenX*cRenderParam.invViewProj[0][2] + screenY*cRenderParam.invViewProj[1][2] + screenZ*cRenderParam.invViewProj[2][2] + cRenderParam.invViewProj[3][2]) / worldW;

	theta[0] = xyz[0]-cRenderParam.eye[0];
	theta[1] = xyz[1]-cRenderParam.eye[1];
	theta[2] = xyz[2]-cRenderParam.eye[2];
	normalize(theta);

	xyz[0] = (xyz[0]-patchXYZ[0])/asi.boxSize[0];
	xyz[1] = (xyz[1]-patchXYZ[1])/asi.boxSize[1];
	xyz[2] = (xyz[2]-patchXYZ[2])/asi.boxSize[2];

	float deltaPoint = 0.0f;
	float errF = 0.0f;


	int its=0;
	if(useRegulaFalsi)
	{
		float b0=cRenderParam.deltaEps*100, a0=-b0; 
		float F, G, w1, fw1;
		EvalBezierBlossom(xyz[0]+a0*theta[0], xyz[1]+a0*theta[1], xyz[2]+a0*theta[2], asi.deg[0], asi.deg[1], asi.deg[2], func, &F);
		EvalBezierBlossom(xyz[0]+b0*theta[0], xyz[1]+b0*theta[1], xyz[2]+b0*theta[2], asi.deg[0], asi.deg[1], asi.deg[2], func, &G);
		its+=2;
		//if(F*G>0)	return;
		while (its<cRenderParam.nIter)
		{
			w1 = (G*a0-F*b0)/(G-F);
			EvalBezierBlossom(xyz[0]+w1*theta[0], xyz[1]+w1*theta[1], xyz[2]+w1*theta[2], asi.deg[0], asi.deg[1], asi.deg[2], func, &fw1);
			if(fabs(fw1)<ERROR_FUNC_EPS)	{xyz[0]+=w1*theta[0], xyz[1]+=w1*theta[1], xyz[2]+=w1*theta[2]; break;}
			if(F*fw1<0)
			{
				b0=w1;
				G=fw1;
				F=F/2.0;
			}
			else
			{
				a0=w1;
				F=fw1;
				G=G/2;
			}
			its++;
		}
		EvalBezierBlossom(xyz[0]+DERIVATIVE_DELTA, xyz[1], xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], func, &nablaF[0]);
		EvalBezierBlossom(xyz[0], xyz[1]+DERIVATIVE_DELTA, xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], func, &nablaF[1]);
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2]+DERIVATIVE_DELTA, asi.deg[0], asi.deg[1], asi.deg[2], func, &nablaF[2]);

		norm[0] = (nablaF[0]-fw1)*DERIVATIVE_DELTA_INV;
		norm[1] = (nablaF[1]-fw1)*DERIVATIVE_DELTA_INV;
		norm[2] = (nablaF[2]-fw1)*DERIVATIVE_DELTA_INV;	
	}
	else
	{
		for(its=1; its<cRenderParam.nIter; its++)
		{
			xyz[0] -= fValue[0];
			xyz[1] -= fValue[1];
			xyz[2] -= fValue[2];

			EvalBezierBlossom(xyz[0], xyz[1], xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], func, &fValue[0]); 
			EvalBezierBlossom(xyz[0]+DERIVATIVE_DELTA, xyz[1], xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], func, &nablaF[0]);
			EvalBezierBlossom(xyz[0], xyz[1]+DERIVATIVE_DELTA, xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], func, &nablaF[1]);
			EvalBezierBlossom(xyz[0], xyz[1], xyz[2]+DERIVATIVE_DELTA, asi.deg[0], asi.deg[1], asi.deg[2], func, &nablaF[2]);

			nablaF[0] = (nablaF[0]-fValue[0])*DERIVATIVE_DELTA_INV;
			nablaF[1] = (nablaF[1]-fValue[0])*DERIVATIVE_DELTA_INV;
			nablaF[2] = (nablaF[2]-fValue[0])*DERIVATIVE_DELTA_INV;


			norm[0] = nablaF[0];
			norm[1] = nablaF[1];
			norm[2] = nablaF[2];

			errF = fValue[0]*fValue[0];

			if ( errF <= ERROR_FUNC_EPS)		break;
			fValue[0] = (fValue[0])/(nablaF[0]*theta[0]+nablaF[1]*theta[1]+nablaF[2]*theta[2]);
			fValue[1] = theta[1]*fValue[0];
			fValue[2] = theta[2]*fValue[0];
			fValue[0] = theta[0]*fValue[0];

			deltaPoint = fabs(fValue[0])+fabs(fValue[1])+fabs(fValue[2]);
			if (deltaPoint <= IntervalEPS)			break;	
		}
		its*=4;
	}

	xyz[0] = xyz[0]*asi.boxSize[0]+patchXYZ[0];
	xyz[1] = xyz[1]*asi.boxSize[1]+patchXYZ[1];
	xyz[2] = xyz[2]*asi.boxSize[2]+patchXYZ[2];

	//float t = (its+0.0)/cRenderParam.nIter;
	//xyz[0] = (1-t)*cRenderParam.material.d[3]+t*cRenderParam.material.d[12];
	//xyz[1] = (1-t)*cRenderParam.material.d[4]+t*cRenderParam.material.d[13];
	//xyz[2] = (1-t)*cRenderParam.material.d[5]+t*cRenderParam.material.d[14];



	OUTPUTVEC(xyz);
	pixel[3] = 1.0f;

	norm[0] /= asi.boxSize[0];
	norm[1] /= asi.boxSize[1];
	norm[2] /= asi.boxSize[2];
	normalize(norm);

	pixel = (float*)(surfaceIn + pixelY*pitch) + 4*pixelX;
	OUTPUTVEC(norm);
}


__global__ void SilAntiAlias(uint* PointOccupied, uint* PointOccupiedValid, uint patchBits,
											unsigned char*surfaceIn, unsigned char*surfaceOut, size_t pitch,
											uint pixelBits,	uint* SilSource, uint actualPoint, float* FuncData, int* PatchIndex, int* PolarIndex,
											ABSSURF_SHARE_INFO asi, int degSize)
{
	uint silPointN = threadIdx.x + blockIdx.x*blockDim.x;
#if MY_DEBUG_CUDA
#endif
	if(silPointN>=actualPoint)							return;

	uint pixelInfo = SilSource[silPointN];
	uint upixelY = (pixelInfo<<(sizeof(uint)*8-pixelBits)>>(sizeof(uint)*8-pixelBits));
	uint upixelX = (pixelInfo>>pixelBits);
	int pixelY = upixelY;
	int pixelX = upixelX;
	float* pixelIn = ((float*)(surfaceIn+pixelY*pitch) + pixelX*4 +2);

	float screenX = 2.0f*(pixelX+0.5f)/cRenderParam.renderWidth-1.0f;
	float screenY = 1.0f-2.0f*(pixelY+0.5f)/cRenderParam.renderHeight;
	float screenZ = *(pixelIn);


	float xyz[3] = {0.0f, 0.0f, 0.0f};
	float norm[3] = {0.0f, 0.0f, 0.0f};
	float worldW;
	int patchIndexXYZ[3];
	float patchXYZ[3];

	float errF = 1000.0f;
	float deltaPoint = 1000.0f;
	float fValue[3] = {0.0f, 0.0f, 0.0f};
	float jacobi[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
	//float norm[3] = {0.0f, 0.0f, 0.0f};
	//int step[3] = {0, 0, 0};
	
	float plane[4] = {0.0f, 0.0f, 0.0f, 0.0f};
	float xyzBak[3];
	float normBak[3];

	worldW = screenX*cRenderParam.invViewProj[0][3] + screenY*cRenderParam.invViewProj[1][3] + screenZ*cRenderParam.invViewProj[2][3] + cRenderParam.invViewProj[3][3];
	xyzBak[0] = xyz[0] = (screenX*cRenderParam.invViewProj[0][0] + screenY*cRenderParam.invViewProj[1][0] + screenZ*cRenderParam.invViewProj[2][0] + cRenderParam.invViewProj[3][0]) / worldW;
	xyzBak[1] = xyz[1] = (screenX*cRenderParam.invViewProj[0][1] + screenY*cRenderParam.invViewProj[1][1] + screenZ*cRenderParam.invViewProj[2][1] + cRenderParam.invViewProj[3][1]) / worldW;
	xyzBak[2] = xyz[2] = (screenX*cRenderParam.invViewProj[0][2] + screenY*cRenderParam.invViewProj[1][2] + screenZ*cRenderParam.invViewProj[2][2] + cRenderParam.invViewProj[3][2]) / worldW;
	
	patchIndexXYZ[0] = floorf((xyz[0]-asi.box0[0])/asi.boxSize[0]);
	patchIndexXYZ[1] = floorf((xyz[1]-asi.box0[1])/asi.boxSize[1]);
	patchIndexXYZ[2] = floorf((xyz[2]-asi.box0[2])/asi.boxSize[2]);

	patchXYZ[0] = patchIndexXYZ[0]*asi.boxSize[0] + asi.box0[0];
	patchXYZ[1] = patchIndexXYZ[1]*asi.boxSize[1] + asi.box0[1];
	patchXYZ[2] = patchIndexXYZ[2]*asi.boxSize[2] + asi.box0[2];

	xyz[0] = (xyz[0]-patchXYZ[0])/asi.boxSize[0];
	xyz[1] = (xyz[1]-patchXYZ[1])/asi.boxSize[1];
	xyz[2] = (xyz[2]-patchXYZ[2])/asi.boxSize[2];


	float EdgeMark = *(pixelIn-2);
	float* pixel = ((float*)(surfaceOut+pixelY*pitch)+ 4*pixelX);
		

	int patchIndex, segInd;
	segInd = patchIndexXYZ[0]*asi.segShift[0] + patchIndexXYZ[1]*asi.segShift[1] + patchIndexXYZ[2]*asi.segShift[2];
	patchIndex = PatchIndex[segInd];
	if(patchIndex<0||PolarIndex[patchIndex]<0)	return;	
	patchIndex *= degSize;


	bool isHor;
	if(IsHorSil(EdgeMark)||IsHorVerSil(EdgeMark))
		isHor = true;
	else if(IsVerSil(EdgeMark))
		isHor = false;
	else return;

	if(isHor)
	{		
		plane[0] = (screenY*cRenderParam.viewProj[0][3]-cRenderParam.viewProj[0][1]);
		plane[1] = (screenY*cRenderParam.viewProj[1][3]-cRenderParam.viewProj[1][1]);
		plane[2] = (screenY*cRenderParam.viewProj[2][3]-cRenderParam.viewProj[2][1]);
		plane[3] = (screenY*cRenderParam.viewProj[3][3]-cRenderParam.viewProj[3][1]) + patchXYZ[0]*plane[0] + patchXYZ[1]*plane[1] + patchXYZ[2]*plane[2];
	}
	else 
	{
		plane[0] = (screenX*cRenderParam.viewProj[0][3]-cRenderParam.viewProj[0][0]);
		plane[1] = (screenX*cRenderParam.viewProj[1][3]-cRenderParam.viewProj[1][0]);
		plane[2] = (screenX*cRenderParam.viewProj[2][3]-cRenderParam.viewProj[2][0]);
		plane[3] = (screenX*cRenderParam.viewProj[3][3]-cRenderParam.viewProj[3][0]) + patchXYZ[0]*plane[0] + patchXYZ[1]*plane[1] + patchXYZ[2]*plane[2];
	}
	
	plane[0] *= asi.boxSize[0];
	plane[1] *= asi.boxSize[1];
	plane[2] *= asi.boxSize[2];


	//int its;
	for(int its=0; its<4; its++)
	{
		xyz[0] -= fValue[0];
		xyz[1] -= fValue[1];
		xyz[2] -= fValue[2];
		
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], FuncData+patchIndex, &fValue[0]);
		EvalBezierBlossom(xyz[0]+DERIVATIVE_DELTA, xyz[1], xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], FuncData+patchIndex, &jacobi[0][0]);
		EvalBezierBlossom(xyz[0], xyz[1]+DERIVATIVE_DELTA, xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], FuncData+patchIndex, &jacobi[0][1]);
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2]+DERIVATIVE_DELTA, asi.deg[0], asi.deg[1], asi.deg[2], FuncData+patchIndex, &jacobi[0][2]);
		jacobi[0][0] = (jacobi[0][0]-fValue[0])*DERIVATIVE_DELTA_INV;
		jacobi[0][1] = (jacobi[0][1]-fValue[0])*DERIVATIVE_DELTA_INV;
		jacobi[0][2] = (jacobi[0][2]-fValue[0])*DERIVATIVE_DELTA_INV;

		norm[0] = jacobi[0][0];
		norm[1] = jacobi[0][1];
		norm[2] = jacobi[0][2];
		if(its==0)
		{
			normBak[0] = norm[0];
			normBak[1] = norm[1];
			normBak[2] = norm[2];
		}

		//EvalBezierBlossomTex(xyz[0], xyz[1], xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], polarTex, patchIndex, &fValue[1]);
		//EvalBezierBlossomTex(xyz[0]+DERIVATIVE_DELTA, xyz[1], xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], polarTex, patchIndex, &jacobi[1][0]);
		//EvalBezierBlossomTex(xyz[0], xyz[1]+DERIVATIVE_DELTA, xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], polarTex, patchIndex, &jacobi[1][1]);
		//EvalBezierBlossomTex(xyz[0], xyz[1], xyz[2]+DERIVATIVE_DELTA, asi.deg[0], asi.deg[1], asi.deg[2], polarTex, patchIndex, &jacobi[1][2]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1], xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], patchIndex, &fValue[1]);
		EvalBezierBlossomPolarTex(xyz[0]+DERIVATIVE_DELTA, xyz[1], xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], patchIndex, &jacobi[1][0]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1]+DERIVATIVE_DELTA, xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], patchIndex, &jacobi[1][1]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1], xyz[2]+DERIVATIVE_DELTA, asi.deg[0], asi.deg[1], asi.deg[2], patchIndex, &jacobi[1][2]);
		jacobi[1][0] = (jacobi[1][0]-fValue[1])*DERIVATIVE_DELTA_INV;
		jacobi[1][1] = (jacobi[1][1]-fValue[1])*DERIVATIVE_DELTA_INV;
		jacobi[1][2] = (jacobi[1][2]-fValue[1])*DERIVATIVE_DELTA_INV;

		fValue[2] = plane[0]*xyz[0] + plane[1]*xyz[1] + plane[2]*xyz[2] + plane[3];

		jacobi[2][0] = plane[0];
		jacobi[2][1] = plane[1];
		jacobi[2][2] = plane[2];

		errF = fValue[0]*fValue[0]+fValue[1]*fValue[1]+fValue[2]*fValue[2];

		if (errF <= ERROR_FUNC_EPS)					break;

		GaussEliminationGPU(jacobi, fValue);


		deltaPoint = fabs(fValue[0])+fabs(fValue[1])+fabs(fValue[2]);
		if (deltaPoint <= IntervalEPS)			break;	
	}

	xyz[0] = xyz[0]*asi.boxSize[0]+patchXYZ[0];
	xyz[1] = xyz[1]*asi.boxSize[1]+patchXYZ[1];
	xyz[2] = xyz[2]*asi.boxSize[2]+patchXYZ[2];
	worldW	   = xyz[0]*cRenderParam.viewProj[0][3] + xyz[1]*cRenderParam.viewProj[1][3] + xyz[2]*cRenderParam.viewProj[2][3] + cRenderParam.viewProj[3][3];
	fValue[0] = (xyz[0]*cRenderParam.viewProj[0][0] + xyz[1]*cRenderParam.viewProj[1][0] + xyz[2]*cRenderParam.viewProj[2][0] + cRenderParam.viewProj[3][0])/worldW; 
	fValue[1] = (xyz[0]*cRenderParam.viewProj[0][1] + xyz[1]*cRenderParam.viewProj[1][1] + xyz[2]*cRenderParam.viewProj[2][1] + cRenderParam.viewProj[3][1])/worldW; 
	fValue[2] = (xyz[0]*cRenderParam.viewProj[0][2] + xyz[1]*cRenderParam.viewProj[1][2] + xyz[2]*cRenderParam.viewProj[2][2] + cRenderParam.viewProj[3][2])/worldW; 

	fValue[0] = ((1.0f+fValue[0])*cRenderParam.renderWidth )/2.0f-0.5f;
	fValue[1] = ((1.0f-fValue[1])*cRenderParam.renderHeight)/2.0f-0.5f;
	//float pixelRefine;
	//if(isHor)		pixelRefine = (fValue[0]-pixelX);
	//else			pixelRefine = (fValue[1]-pixelY);

	float pixelNew;
	int pixelOri;
	//int pixelRefine, pixelIter, dir, lim, dir4;
	if(isHor)
	{		
		pixelNew = fValue[0];
		pixelOri = pixelX;
		//dir = (pixelNew>pixelOri) ? 1 : -1;
		//dir4 = 4*dir;
		//lim = cRenderParam.renderWidth;
	}
	else
	{
		pixelNew = fValue[1];
		pixelOri = pixelY;
		//dir = (pixelNew>pixelOri) ? 1 : -1;
		//dir4 = int(pitch/4)*dir;
		//lim = cRenderParam.renderHeight;
	}

	float* pixelNorm = pixelIn-2;
	//if(errF<=ERROR_FUNC_EPS&&fabs(pixelNew-pixelOri)<3/*&&fabs((screenZ-fValue[2])/screenZ)<0.01*/)
	//{
	//	norm[0] /= asi.boxSize[0];
	//	norm[1] /= asi.boxSize[1];
	//	norm[2] /= asi.boxSize[2];
	//	normalize(norm);

	//	float pad[3], padN[3], stepN;

		//if(dir>0)	pixelRefine = int(ceil(pixelNew));
		//else		pixelRefine = int(floor(pixelNew));
		//pixelRefine = min(pixelRefine, lim-1);
		//pixelRefine = max(0, pixelRefine);

	//	pixel-=dir4;
	//	pixelNorm-=dir4;
	//	if((pixelX-dir)>=0&&(pixelX-dir)<lim)
	//	{
	//		stepN = 1.0/(abs(pixelRefine-pixelOri)+1);
	//		xyzBak[0] = pixel[0];
	//		xyzBak[1] = pixel[1];
	//		xyzBak[2] = pixel[2];
	//		pad[0] = (xyz[0]-xyzBak[0])*stepN;
	//		pad[1] = (xyz[1]-xyzBak[1])*stepN;
	//		pad[2] = (xyz[2]-xyzBak[2])*stepN;
	//		normBak[0] = pixelNorm[0];
	//		normBak[1] = pixelNorm[1];
	//		normBak[2] = pixelNorm[2];
	//		padN[0] = (norm[0]-normBak[0])*stepN;
	//		padN[1] = (norm[1]-normBak[1])*stepN;
	//		padN[2] = (norm[2]-normBak[2])*stepN;
	//	}
	//	else	return;		

	//	for(pixelIter=pixelOri; (pixelIter!=pixelRefine); pixelIter+=dir, pixel+=dir4, pixelNorm+=dir4)
	//	{
	//		xyzBak[0] += pad[0];
	//		xyzBak[1] += pad[1];
	//		xyzBak[2] += pad[2];
	//		pixel[0] = xyzBak[0];
	//		pixel[1] = xyzBak[1];
	//		pixel[2] = xyzBak[2];
	//		pixel[3] = 1.0f;
	//		
	//		normBak[0] += padN[0];
	//		normBak[1] += padN[1];
	//		normBak[2] += padN[2];
	//		pixelNorm[0] = normBak[0];
	//		pixelNorm[1] = normBak[1];
	//		pixelNorm[2] = normBak[2];
	//	}
	//
	//	OUTPUTVEC(xyz);
	//	pixel[3] = (pixelNew-pixelRefine)*dir;

	//
	//	pixelNorm[0] = norm[0];
	//	pixelNorm[1] = norm[1];
	//	pixelNorm[2] = norm[2];
	//	pixelNorm[3] = isHor;
	//}
	//else
	{
		OUTPUTVEC(xyzBak);
		pixel[3] = 0.5/*(pixelNew-pixelOri)*/;

		normBak[0] /= asi.boxSize[0];
		normBak[1] /= asi.boxSize[1];
		normBak[2] /= asi.boxSize[2];
		normalize(normBak);
		pixelNorm[0] = normBak[0];
		pixelNorm[1] = normBak[1];
		pixelNorm[2] = normBak[2];
		pixelNorm[3] = (int)isHor;
	}
}

extern "C"
void MultiSilRender0726(RENDER_PARAM* multiSilRender, void* surfaceIn, void* surfaceOut, size_t pitch, int iter)
{
	cudaMemcpyToSymbol(cRenderParam, multiSilRender, sizeof(RENDER_PARAM));

	CUDA_SAFE_CALL( cudaMemset(dPolarIndex, -1, nPatch*sizeof(int)) );
	float3 eye = make_float3(multiSilRender->eye[0], multiSilRender->eye[1], multiSilRender->eye[2]);
	//int degSize = (asi.deg[0]+1)*(asi.deg[1]+1)*(asi.deg[2]+1); 
	size_t Ns = degSize*20*sizeof(float);

	GetPolarData<<<nPatch, degSize, Ns>>>(dFuncData, dPolarData, /*dPatchIndex, */dPatchIndexInverse, dPolarIndex, eye, asi);
	//dumpBufferInt(dPolarIndex, nPatch);


#if DEBUG_BUFFERS
	printf("dPolarData:\n");
	dumpBufferfloat(dPolarData, degSize*nPatch);
#endif

	//int silSize = renderWidth*renderHeight*sizeof(uint);
	CUDA_SAFE_CALL( cudaMemset(dSilOccupied, 0, allPointSize) );
	CUDA_SAFE_CALL( cudaMemset(dSilOccupiedValid, 0, allPointSize) );
	dim3 Db = dim3(16,16);
	dim3 Dg = dim3( renderWidth/Db.x+1, renderHeight/Db.y+1 );

	CUDA_SAFE_CALL( cudaMemset(dPointOccupied, 0, allPointSize) );
	CUDA_SAFE_CALL( cudaMemset(dPointOccupiedValid, 0, allPointSize) );
	CheckAllSilPointUniform<<<Dg,Db>>>((unsigned char*)surfaceIn, pitch, pixelBits, dSilOccupied, dSilOccupiedValid, patchBits,
								dPatchIndex, dPointOccupied, dPointOccupiedValid, asi);


#if DEBUG_BUFFERS
	printf("dSilOccupied:\n");
	dumpBuffer(dPointOccupiedValid, renderWidth*renderHeight);
#endif

	uint compNum;

#if DEBUG_BUFFERS
	printf("Debug Information:\n");
	dumpBuffer((uint*)dDebugInfo, silNum);
#endif

	uint sortPadNum;
	compact(&compNum, allPointSize, dNumber, dPointOccupied, dPointOccupiedValid, dOccupiedCompact, 0, NULL);   
	sortPadAll(compNum, &sortPadNum, allPointSize, dOccupiedCompact, dOccupiedSort, dOccupiedSortPad);	
	Db = dim3(PIXELTHREADS, 1, 1);
	Dg = dim3(sortPadNum/PIXELTHREADS, 1, 1);
#ifdef NORMALFUNCTION
	Ns = PIXELTHREADS*sizeof(float)*2;
#else
	//512&1024 work fine with deg.xyz=3;
	Ns =1024; //(asi.deg[0]+1)*degShift[0]*3*sizeof(float)*2
#endif
	//PointOccupied for error, SilOccupied for Boundary.
	CUDA_SAFE_CALL( cudaMemset(dPointOccupied, 0, allPointSize) );
	CUDA_SAFE_CALL( cudaMemset(dPointOccupiedValid, 0, allPointSize) );

	Newton<<<Dg, Db, Ns>>>((unsigned char*)surfaceIn, (unsigned char*)surfaceOut, pitch, pixelBits, patchBits, 
		dOccupiedSortPad, dPointOccupied, dPointOccupiedValid, dSilOccupied, dSilOccupiedValid, dFuncData, dPatchIndex, dPatchIndexInverse, dPolarIndex, asi, iter);

	uint silNum;
	//dumpBuffer(dSilOccupied, renderWidth*renderHeight);

	compact(&silNum, allPointSize, dNumber, dSilOccupied, dSilOccupiedValid, dOccupiedCompact, 0, NULL/*renderWidth*15, dSil*/);
	//printf("%u sil points.\n", silNum);


	SilAntiAlias<<<silNum/PIXELTHREADS+1, PIXELTHREADS>>>(dPointOccupied, dPointOccupiedValid, patchBits,
								(unsigned char*)surfaceIn, (unsigned char*)surfaceOut, pitch, pixelBits,
								dOccupiedCompact, silNum, dFuncData, dPatchIndex, dPolarIndex, asi, degSize);	

#if DEBUG_BUFFERS
	printf("dPointOccupiedValid:\n");
	dumpBuffer(dOccupiedCompact, maxPoint);
#endif
}